@jorgepiloto had suggested plotting the satellite visibility function values to check the behavior at the two suspected (lat, on) values.

For the case of (46 deg, 45 deg) for (lat, on) values, which gives an error, I got:

(Here, the x-axis is not the time but shows "k", where the __call__ method is called kth time.)

With terminal=True:

With terminal=False

poliastro's event detector detects the event at: 2020-01-01 02:00:37.461 whereas orekit detects it at 2020-01-01 03:33:36.580.

For a (lat, lon) pair which doesn't yield an error, there are only two instances of "weird spikes" instead of the three spikes seen for the mismatch case. I think the time mismatch is because of this.

To clarify, this is the diff:

diff --git a/src/poliastro/core/events.py b/src/poliastro/core/events.py
index 483aede3..66d9fca4 100644
--- a/src/poliastro/core/events.py
+++ b/src/poliastro/core/events.py
@@ -102,7 +102,10 @@ def visibility_function(k, u_, phi, theta, R, R_p, H):
         ]
     )
     new_rho = rot_matrix @ rho
+    print("New rho before normalization: ", new_rho)
     new_rho = new_rho / np.linalg.norm(new_rho)
+    print("New rho after normalization: ", new_rho)
     el = np.arcsin(new_rho[-1])
+    print("Elevation: ", el)
 
     return el
diff --git a/src/poliastro/twobody/events.py b/src/poliastro/twobody/events.py
index 3bd57cde..4d8e6278 100644
--- a/src/poliastro/twobody/events.py
+++ b/src/poliastro/twobody/events.py
@@ -279,6 +279,7 @@ class SatelliteVisibilityEvent(Event):
         self._R = orbit.attractor.R.to(u.km).value
         self._k = orbit.attractor.k.to_value(u.km ** 3 / u.s ** 2)
         self._R_p = orbit.attractor.R_polar.to(u.km).value
+        self.l = []
 
     def __call__(self, t, u_, k):
         self._last_t = t
@@ -303,4 +304,5 @@ class SatelliteVisibilityEvent(Event):
             self._R_p,
             self._H.to_value(u.km),
         )
+        self.l.append(el)
         return el

this is a script to reproduce the spikes:

import numpy as np
import pytest
from astropy import units as u
from astropy.tests.helper import assert_quantity_allclose
from astropy.time import Time
from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.twobody.events import (
    LatitudeCrossEvent,
    NodeCrossEvent,
    PenumbraEvent,
    SatelliteVisibilityEvent,
    UmbraEvent,
)
from poliastro.twobody.propagation import cowell

import matplotlib.pyplot as plt

# Macro for conversion between degrees and radians
DEG_TO_RAD = np.pi / 180

a, ecc, inc, raan, argp, nu = 6828137.0, 0.4, 87.0, 20.0, 10.0, 0

epoch0_poliastro = Time("2020-01-01", scale="utc")
ss0_poliastro = Orbit.from_classical(
    Earth,
    a * u.m,
    ecc * u.one,
    inc * u.deg,
    raan * u.deg,
    argp * u.deg,
    nu * u.deg,
    epoch0_poliastro,
)

# Unpack orekit and poliastro events
poliastro_event = SatelliteVisibilityEvent(
    ss0_poliastro, 46 * u.deg, 45 * u.deg, 5.0 * u.m, terminal=False,
)

# Time of fliht for propagating the orbit
tof = (1 * u.day).to_value(u.s)

# Propagate poliastro's orbit
_, _ = cowell(
    Earth.k,
    ss0_poliastro.r,
    ss0_poliastro.v,
    # Generate a set of tofs, each one for each propagation second
    np.linspace(0, tof, 100) * u.s,
    events=[poliastro_event],
)
poliastro_event_epoch = ss0_poliastro.epoch + poliastro_event.last_t

# Test both event epochs by checking the distance in seconds between them
fig, ax = plt.subplots()
ax.plot(poliastro_event.l, "x")
ax.axhline(linestyle="--", color="k")
ax.set_ylabel("Visibility function value")
ax.set_title("Satellite visibility (terminal=False)")
fig.savefig("vis.png")

and this is how the plot looks like:

The spikes are to be expected, they are part of SciPy trying to find the crossing point of the event.

Hey @Yash-10, can you plot the discrete elevation values around the time where you are experiencing trouble? It would be good to compare them for orekit and poliastro. I suspect you are hitting a borderline case with very low elevation and the issue probably does not help.

The one time huge (5500 sec) error suggests me that you are off by a full orbital period, no? (didn't calculate it, but looks a bit like that to me)

The one time huge (5500 sec) error suggests me that you are off by a full orbital period, no? (didn't calculate it, but looks a bit like that to me)

Thanks, @egemenimre for this point, and sorry for the delay. This is very interesting since indeed it almost seems to happen. I observed the behavior again and saw the following:

• The time period of the orbit is 5615.18823984 s and the discrepancy is 5289.606 s.
• poliastro's event is detected after 1333.82 s whereas orekit's is detected after 6623.426 s. So poliastro detected the event around one orbital period before the actual occurrence of the event.

Since these elevation values are calculated at time instants decided by SciPy, I think it is a bit involved to plot elevation as a function of time. I would try to do it and update it here. The above satellite visibility function plot (with the y-axis showing the elevation, although the x-axis not showing the time but the count of time instances (decided internally by SciPy)) seems to have already been created by @astrojuanlu above.

Hi there. Actually time is not strictly required for the elevation plot. What I am interested in is to see how close the elevation is to zero around this peculiar incident. My hypothesis is that, at that specific time there is a near-horizon event (let's say the satellite climbs only half a degree over the horizon and then dips). As the event finder algorithm probably looks for the near-zero value, and does some weird stuff trying to zero in. In a high elevation pass, the zero crossing is very clean but with a very low elevation pass with discrete points, finding the zero crossing can be very tricky (and inaccurate). I suspect SciPy might be struggling there.

The other point is that, this is a bit of a dangerous method to check for the pass times. Let's say Orekit finds 10 events. But because of a low elevation pass, poliastro decides that there has been no pass for pass number 3. As this call subtracts the events one by one, orekit event 3 will correspond to poliastro event 3 (which is actually orekit event 4). And then you will have this mismatch until the end.

# Test both event epochs by checking the distance in seconds between them
dt = np.abs((orekit_event_epoch - poliastro_event_epoch).to(u.s))

Oh, now I understand. It looks like SciPy detects an event if the visibility function is very close but perturbing around zero whereas orekit doesn't consider such small perturbations around zero but only the larger perturbations. The below plots could make that clearer, I think:

A) A case where poliastro and orekit do not match

i) terminal=False

Notice the region (magnified) below from the above plot:

It feels SciPy is trying hard in that region.

ii) terminal=True

B) A case where poliastro and orekit match

i) terminal=False

ii) terminal=True

I will take a step back and point out to an issue (to @astrojuanlu as well). Admittedly I have worked on event detection stuff post-facto, meaning that you run a propagation for the full duration and then use the trajectory to find the "communications opportunities" events. For GS pass calculation this usecase is very common. The way Orekit (and now poliastro) handles things through the internal event handling mechanism of the ODE Solver may be more suited to tasks in trajectory computation and optimisation ("when altitude drops to 400 km, compute a firing" or "stop electric propulsion when in eclipse"). I'm not saying the approach is in any way wrong, all I'm saying is that pass computation may be a relatively less common usecase for this function in reality. The post facto approach means that I have the full trajectory and I know many points around the low elevation pass. To this I can fit a polynomial (or more accurately, a bunch of splines) and I can compute the "roots" of that polynomial (zero crossings) rather easily.

Coming back to the original point, I'm not sure about what you mean by perturbations - do you mean the back and forth of the internal event handler to find the zero crossing? I suspect the internal event handler does something like this:

1. the value is negative, all good, let's continue the propagation
2. ah the value turned positive in the next step, there must be a zero crossing event here somewhere
3. let me go back and compute things in a smaller stepsize until I converge on the exact spot where the zero-crossing happened

Your plots show clearly Scipy is hopping above and below the zero line, to find the exact crossing time. Case (B ii) is showing the 3 steps done right. Case (A i) shows the struggle of the event finder going up and down to find the zero crossing.

Now, if your satellite has a high elevation pass, the operations above are clean and accurate. However, the low elevation pass problem is a very difficult problem to solve like that. You are trying to find the existence and then the begin and end times of a tiny part of parabola.

Due to numerical issues, the algorithm can easily get lost within the 3 steps above and decide that, actually, there was no pass because there has been no zero crossing (as a crude example, you could see this with very large timesteps, you would simply jump over the very brief elevation >0 event). Similarly, the ODE solver will suffer a bit when you have an event where the elevation comes very close to zero, but not exceeds it.

This problem is similar to finding the exact time in the day where the Sun is at its highest. You look for the derivative = 0, but it is changing so slowly, it is not very easy to calculate it very accurately.

I had a similar problem in my (post-facto) eclipse finder and I found the solution in changing the equivalent of your visibility function.

So far my assumption has been Orekit and poliastro propagations yield exactly the same results and coordinate frame conversions yield exactly the same results. This should also be tested. Any difference in propagation (force model as well as the numerical model) may cause a very low elevation pass event to disappear in the other. Similarly, if the coordinate conversions use different IERS EOP files, then the results would also be different. So I would make sure there aren't any issues with these two. This is why I would simply check the difference between Poliastro and Orekit in the following:

1. Ground location coords in ITRF
2. Satellite coords in ITRF
3. Elevation

And maybe these questions might help:

1. Do you see any difference in the event times when you have terminal on or off?
2. Do you see any difference in event times when you manipulate stepsize through tolerance (atol and rtol)?

The final and the most unfortunate conclusion might be that, assuming Orekit results are the truth, the Hipparchus ODE event handling mechanism may be simply better than the one in Scipy.

It had been a long post, but this is a difficult problem. The event handler in the ODE Solver is opaque and you are trying to compare two completely different systems (Orekit and Astropy/Scipy based poliastro) where you can have a lot of inconsistencies.

It makes more sense now what the underlying issue might be.

I'm not sure about what you mean by perturbations - do you mean the back and forth of the internal event handler to find the zero crossing?

Yeah, that's what I meant, sorry for using "perturbation" vaguely (especially in orbital mechanics!).

Your three-step explanation of the underlying event solver made me curious to experiment with the different solvers already available in SciPy (i.e. the method argument to solve_ivp). The default poliastro uses is DOP853 and the documentation states:

‘DOP853’ is recommended for solving with high precision (low values of rtol and atol).

I think this might be the root cause of the observed mismatches, since by using DOP853, poliastro's implementation is being sensitive with low elevation values, as you state.

Interesting observation

This stems from your point:

Do you see any difference in event times when you manipulate stepsize through tolerance (atol and rtol)?

To prevent poliastro to detect an event at the low-elevation edge case, I tried to do the following:

• Use the LSODA method (used LSODA since SciPy could then allow adding a min_step argument to solve_ivp which would not be otherwise possible for other solver methods).
• Add a min_step argument to solve_ivp by setting it to 1e-1, for example. This would ensure no step size is smaller than 1e-1.

What I observed after doing this and running the validation again, the mismatch reduces from the earlier ~5500 s to ~2 min

Although this might seem a good direction to solve such a problem, by manipulating the step size, we might hamper some/all solutions that were agreeing well with orekit earlier.

Maybe implementing a custom ODESolver and using it inside solve_ivp instead could allow more control, but I note the caveats about ODE event handlers you mentioned.

Small updates:

• I confirm that the added implementation matches with example 5.9 in terms of final elevation values and all other intermediate results.
• Doing

print(ss0_orekit.getPVCoordinates(topo_frame))
print(ss0_poliastro.r, ss0_poliastro.v)

yields different values whereas doing

print(ss0_orekit.getPVCoordinates())
print(ss0_poliastro.r, ss0_poliastro.v)

yields the same values.

• Satellite coordinates expected by the implementation are geocentric equatorial. Question is which system does poliastro and orekit assume for the position vectors?

I am taking a look at this after a few months. I again compared this satellite visibility implementation from Howard Curtis fourth edition example 5.9 (see poliastro/poliastro#1299) and it matches with it. So, the implementation is in-line with that example.

Comparing elevations from poliastro and orekit by plotting on the same graph seems a very difficult thing to do from a user-level since both event detection mechanisms work in a very different way. However, as suggested by @egemenimre, I have tried to do some visualizations that could be done by manipulating the script shared by @astrojuanlu above.

(This script requires orekit and could be run from the directory: validation/orekit/tests. It involves some hackery for plotting the values.

Disclaimer: The graphs should not be directly compared since the corresponding values are not taken at the same corresponding time instants and these plots only give a rough view of how the elevation value changes.)

import numpy as np
import pytest
from astropy import units as u
from astropy.tests.helper import assert_quantity_allclose
from astropy.time import Time
from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.twobody.events import (
    LatitudeCrossEvent,
    NodeCrossEvent,
    PenumbraEvent,
    SatelliteVisibilityEvent,
    UmbraEvent,
)
from poliastro.twobody.propagation import cowell

import matplotlib.pyplot as plt

# Orekit imports
import orekit
from orekit.pyhelpers import setup_orekit_curdir

# Setup orekit virtual machine and associated data
VM = orekit.initVM()
setup_orekit_curdir("orekit-data.zip")

from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.propagation.analytical import KeplerianPropagator
from org.orekit.propagation import SpacecraftState
from org.orekit.orbits import KeplerianOrbit, PositionAngle
from org.orekit.utils import Constants, IERSConventions
from org.orekit.frames import FramesFactory, TopocentricFrame
from org.orekit.bodies import GeodeticPoint, OneAxisEllipsoid
from org.orekit.propagation.events import ElevationDetector
from org.orekit.propagation.events.handlers import (
    # StopOnDecreasing,
    StopOnEvent,
    # StopOnIncreasing,
)

# Macro for conversion between degrees and radians
DEG_TO_RAD = np.pi / 180

a, ecc, inc, raan, argp, nu = 6828137.0, 0.0073, 87.0, 20.0, 10.0, 0

epoch0_poliastro = Time("2020-01-01", scale="utc")
ss0_poliastro = Orbit.from_classical(
    Earth,
    a * u.m,
    ecc * u.one,
    inc * u.deg,
    raan * u.deg,
    argp * u.deg,
    nu * u.deg,
    epoch0_poliastro,
)

# Unpack orekit and poliastro events
poliastro_event = SatelliteVisibilityEvent(
    ss0_poliastro, 45 * u.deg, 47 * u.deg, 5.0 * u.m, terminal=True,
)

# Time of fliht for propagating the orbit
tof = float((1 * u.day).to_value(u.s))

# Propagate poliastro's orbit
_, _ = cowell(
    Earth.k,
    ss0_poliastro.r,
    ss0_poliastro.v,
    # Generate a set of tofs, each one for each propagation second
    np.linspace(0, tof, 100) * u.s,
    events=[poliastro_event],
)
poliastro_event_epoch = ss0_poliastro.epoch + poliastro_event.last_t

### Orekit ###
epoch0_orekit = AbsoluteDate(2020, 1, 1, 0, 0, 00.000, TimeScalesFactory.getUTC())
ss0_orekit = KeplerianOrbit(
    float(a),
    float(ecc),
    float(inc * DEG_TO_RAD),
    float(argp * DEG_TO_RAD),
    float(raan * DEG_TO_RAD),
    float(nu * DEG_TO_RAD),
    PositionAngle.TRUE,
    FramesFactory.getEME2000(),
    epoch0_orekit,
    Constants.WGS84_EARTH_MU,
)
# Define a flattened Earth instead of a perfect sphere model for the satellite visibility event.
Earth_orekit_satellite_visibility = OneAxisEllipsoid(
    Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
    Constants.WGS84_EARTH_FLATTENING,
    FramesFactory.getITRF(IERSConventions.IERS_2010, True),
)
point = GeodeticPoint(45 * DEG_TO_RAD, 47 * DEG_TO_RAD, 5.0)  # lat, lon, h
topo_frame = TopocentricFrame(Earth_orekit_satellite_visibility, point, "earth-station")
orekit_event = ElevationDetector(topo_frame).withHandler(StopOnEvent())
orekit_event2 = ElevationDetector(topo_frame).withHandler(StopOnEvent())

# Build the orekit propagator and add the event to it.
propagator = KeplerianPropagator(ss0_orekit)
propagator2 = KeplerianPropagator(ss0_orekit)
propagator.addEventDetector(orekit_event)
propagator2.addEventDetector(orekit_event2)
orekit_elevations = []
for i in range(0, int(tof), 1140):
    state = propagator.propagate(epoch0_orekit.shiftedBy(float(i)), epoch0_orekit.shiftedBy(float(i)).shiftedBy(float(1)))
    state = SpacecraftState(state.orbit)
    elevation = orekit_event.g(state)
    orekit_elevations.append(elevation)

orekit_event_epoch_raw = state.orbit.getPVCoordinates().getDate()
# Convert orekit epoch to astropy Time instance
orekit_event_epoch_str = orekit_event_epoch_raw.toString(TimeScalesFactory.getUTC())
orekit_event_epoch = Time(orekit_event_epoch_str, scale="utc", format="isot")
orekit_event_epoch.format = "iso"

state2 = propagator2.propagate(epoch0_orekit, epoch0_orekit.shiftedBy(tof))
orekit_event_epoch_raw2 = state2.orbit.getPVCoordinates().getDate()
# Convert orekit epoch to astropy Time instance
orekit_event_epoch_str2 = orekit_event_epoch_raw2.toString(TimeScalesFactory.getUTC())
orekit_event_epoch2 = Time(orekit_event_epoch_str2, scale="utc", format="isot")
orekit_event_epoch2.format = "iso"

poliastro_event_epoch = ss0_poliastro.epoch + poliastro_event.last_t

print(poliastro_event_epoch)
print(orekit_event_epoch)

# Test both event epochs by checking the distance in seconds between them
fig, ax = plt.subplots()
ax.plot(poliastro_event.l, label='poliastro')
ax.plot(orekit_elevations, label='orekit')
ax.legend()
ax.axhline(linestyle="--", color="k")
ax.set_ylabel("Visibility function value")
fig.savefig("vis.png")
plt.show()

which gives this plot:

(I belive there is no need to worry about the fact that orekit's plot does not even reach zero till the end because I only sampled the values after every 1 sec. In reality, it would have reached zero)

and the output print is:

poliastro: 2020-01-01 02:00:46.712
orekit: 2020-01-01 23:45:01.000

There's a 21-22 hour mismatch. From the plot, it is clearly seen that poliastro and orekit both reach zero around the x-axis value of ~45. It seems that orekit was just not sensitive to that because of which orekit detects it some later time?